############### Appendix ################
rm(list=ls())
gc()

library(dplyr)
library(rstan)
library(ggpubr)

######### 1. Evaluation of IRT models #############

######### 1.1. Party competition space

### Static model  - Figures A1, A2, and A3
load("model/Model party ideal position.rda")

model_iter <- extract(model1)

## Figure A1

p1.1 <- ggplot()+
  geom_line(aes(y=model_iter$theta[,5], x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

p1.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$theta[1:10000,5], model_iter$theta[10001:15000,5] * -1, model_iter$theta[15001:20000,5]), x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

png("Plots/Figure A1.png", width = 3000, height=1500, res=300)
ggarrange(p1.1, p1.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A2

p2.1 <- ggplot()+
  geom_line(aes(y=model_iter$beta[,429], x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-12,12))+
  theme(legend.position = "bottom")

p2.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$beta[1:10000,429], model_iter$beta[10001:15000,429] * -1, model_iter$beta[15001:20000,429]), x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[140]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-12,12))+
  theme(legend.position = "bottom")



png("Plots/Figure A2.png", width = 3000, height=1500, res=300)
ggarrange(p2.1, p2.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A3


png("Plots/Figure A3.png", width = 3000, height=1500, res=300)
ggplot()+
  geom_line(aes(y=model_iter$alpha[,429], x = rep(5001:10000, 4), color = factor(rep(1:4, 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(a[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")
dev.off()

### Dynamic model - Figures A4, A5, and A6
load("model/Model party ideal position_by legislature.rda")


model_iter <- extract(model2)

# Figure A4

p4.1 <- ggplot()+
  geom_line(aes(y=model_iter$theta2[,16,5], x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta["svp,t=16"]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

p4.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$theta2[1:10000,16,5], model_iter$theta2[10001:15000,16,5] * -1, model_iter$theta2[15001:20000,16,5]), x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

png("Plots/Figure A4.png", width = 3000, height=1500, res=300)
ggarrange(p4.1, p4.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A5

p5.1 <- ggplot()+
  geom_line(aes(y=model_iter$beta[,429], x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-12,12))+
  theme(legend.position = "bottom")

p5.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$beta[1:10000,429], model_iter$beta[10001:15000,429] * -1, model_iter$beta[15001:20000,429]), x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-12,12))+
  theme(legend.position = "bottom")



png("Plots/Figure A5.png", width = 3000, height=1500, res=300)
ggarrange(p5.1, p5.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()


## Figure A6

png("Plots/Figure A6.png", width = 3000, height=1500, res=300)
ggplot()+
  geom_line(aes(y=model_iter$alpha[,429], x = rep(5001:10000, 4), color = factor(rep(1:4, 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(a[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")
dev.off()



######### 1.2. Party canton space - Figures A7, A8, and A9
load("model/Model IRT with parties and cantons.Rda")

thetas_party_cantons <- summary(model_party_canton, pars = "theta")
thetas_party_cantons <- as.data.frame(thetas_party_cantons)

ggplot()+
  geom_point(aes(x=thetas_party_cantons$c_summary.50..chain.1, y=1:31, color = c(rep("black", 26), rep("grey", 5))))


model_iter <- extract(model_party_canton)

## Figure A7

p7.1 <- ggplot()+
  geom_line(aes(y=model_iter$theta[,31], x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

p7.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$theta[1:10000,31], model_iter$theta[10001:15000,31] * -1, model_iter$theta[15001:20000,31]), x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

png("Plots/Figure A7.png", width = 3000, height=1500, res=300)
ggarrange(p7.1, p7.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A8

p8.1 <- ggplot()+
  geom_line(aes(y=model_iter$beta[,429], x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-12,12))+
  theme(legend.position = "bottom")

p8.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$beta[1:10000,429], model_iter$beta[10001:15000,429] * -1, model_iter$beta[15001:20000,429]), x = rep(5001:10000, 4), color = factor(rep(1:4, each = 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-12,12))+
  theme(legend.position = "bottom")



png("Plots/Figure A8.png", width = 3000, height=1500, res=300)
ggarrange(p8.1, p8.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A9

png("Plots/Figure A9.png", width = 3000, height=1500, res=300)
ggplot()+
  geom_line(aes(y=model_iter$alpha[,429], x = rep(5001:10000, 4), color = factor(rep(1:4, 5000))))+
  xlab("Iteration after warmups")+
  ylab(expression(a[140]))+
  labs(color = "Chains")+
  theme_minimal()+
  theme(legend.position = "bottom")
dev.off()



######### 1.3. Party municipality space - Figures A10, A11, and A12
load("model/Model IRT with parties and municipalities.Rda")


model_iter <- extract(model_party_mun)

## Figure A10

p10.1 <- ggplot()+
  geom_line(aes(y=model_iter$theta[,2177], x = rep(501:1000, 4), color = factor(rep(1:4, each = 500))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

p10.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$theta[1:1000,2177], model_iter$theta[1001:2000,2177] * -1), x = rep(501:1000, 4), color = factor(rep(1:4, each = 500))))+
  xlab("Iteration after warmups")+
  ylab(expression(theta[svp]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-5,5))+
  theme(legend.position = "bottom")

png("Plots/Figure A10.png", width = 3000, height=1500, res=300)
ggarrange(p10.1, p10.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A11

p11.1 <- ggplot()+
  geom_line(aes(y=model_iter$beta[,429], x = rep(501:1000, 4), color = factor(rep(1:4, each = 500))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-2,2))+
  theme(legend.position = "bottom")

p11.2 <- ggplot()+
  geom_line(aes(y=c(model_iter$beta[1:1000,429], model_iter$beta[1001:2000,429] * -1), x = rep(501:1000, 4), color = factor(rep(1:4, each = 500))))+
  xlab("Iteration after warmups")+
  ylab(expression(b[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-2,2))+
  theme(legend.position = "bottom")



png("Plots/Figure A11.png", width = 3000, height=1500, res=300)
ggarrange(p11.1, p11.2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()

## Figure A12

png("Plots/Figure A12.png", width = 3000, height=1500, res=300)
ggplot()+
  geom_line(aes(y=model_iter$alpha[,429], x = rep(501:1000, 4), color = factor(rep(1:4, 500))))+
  xlab("Iteration after warmups")+
  ylab(expression(a[429]))+
  labs(color = "Chains")+
  theme_minimal()+
  ylim(c(-2,2))+
  theme(legend.position = "bottom")
dev.off()


############## Static VS dynamic party competition space #################

## Figure A13

load("model/Model party ideal position.rda")
load("model/Model party ideal position_by legislature.rda")


betas_static <- summary(model1, pars = "beta")
betas_dynamic <- summary(model2, pars = "beta")

betas_static <- as.data.frame(betas_static)
betas_dynamic <- as.data.frame(betas_dynamic)

png("Plots/Figure A13.png", width=2000, height=2000, res=300)
ggplot()+
  geom_point(aes(x = betas_static$c_summary.50..chain.1, y = betas_dynamic$c_summary.50..chain.1))+
  geom_errorbarh(aes(xmin = betas_static$c_summary.50..chain.1 - 1.96 * betas_static$c_summary.sd.chain.1, 
                     xmax = betas_static$c_summary.50..chain.1 + 1.96 * betas_static$c_summary.sd.chain.1, 
                     y = betas_dynamic$c_summary.50..chain.1), alpha = .2)+
  geom_errorbar(aes(ymin = betas_dynamic$c_summary.50..chain.1 - 1.96 * betas_dynamic$c_summary.sd.chain.1, 
                     ymax = betas_dynamic$c_summary.50..chain.1 +1.96 * betas_dynamic$c_summary.sd.chain.1, 
                     x = betas_static$c_summary.50..chain.1), alpha = .2)+
  theme_minimal()+
  xlab("Discrimination of ballots in the Static Party Competition Space")+
  ylab("Discrimination of ballots in the Dyanmic Party Competition Space")+
  geom_text(aes(x=8, y=-10, label = paste0("Pearson cor = ", round(cor(betas_static$c_summary.50..chain.1, betas_dynamic$c_summary.50..chain.1), 3))))
dev.off()

## Figure A14


load("data/Data model cantons and municipalities.Rda")

data_canton_scale_dynamic <- transform(data_canton_mun$data_canton_simple, 
                                       Discrimination = scale(c_summary.50..chain.1),
                                       Discrimination.2.5 = scale(c_summary.2.5..chain.1),
                                       Discrimination.25 = scale(c_summary.25..chain.1),
                                       Discrimination.75 = scale(c_summary.75..chain.1),
                                       Discrimination.97.5 = scale(c_summary.97.5..chain.1),
                                       diff_log_sc = scale(diff_log),
                                       log_sc = log_YesShare)


discrimination <- data_canton_scale_dynamic[data_canton_scale_dynamic$Place=="- Zug",]

discrimination <- discrimination %>% 
  group_by(year) %>% 
  summarise(disc = mean(abs(c_summary.50..chain.1)))


png("Plots/Figure A14.png", width=3000, height=2000, res=300)
ggplot(discrimination)+
  geom_point(aes(x=year, y = disc))+
  labs(y="Average ballots' discrimination absolute value", x="Year")+
  theme_minimal()
dev.off()

## Figure A15

load("data/Results position municipalities and cantons.Rda")

data_canton_results <- data_results_canton_mun$data_canton_results

png("Plots/Figure A15.png", width=2500, height=2500, res = 300)
ggarrange(
  ggplot(data_canton_results)+
    geom_point(aes(x=rand, y = rand.2.5))+
    geom_errorbar(aes(x=rand, ymin = low.2.5, ymax = high.2.5))+
    geom_errorbarh(aes(y=rand.2.5, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i2.5])),
  
  ggplot(data_canton_results)+
    geom_point(aes(x=rand, y = rand.25))+
    geom_errorbar(aes(x=rand, ymin = low.25, ymax = high.25))+
    geom_errorbarh(aes(y=rand.25, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i25])),
  
  ggplot(data_canton_results)+
    geom_point(aes(x=rand, y = rand.75))+
    geom_errorbar(aes(x=rand, ymin = low.75, ymax = high.75))+
    geom_errorbarh(aes(y=rand.75, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i75])),
  
  ggplot(data_canton_results)+
    geom_point(aes(x=rand, y = rand.97.5))+
    geom_errorbar(aes(x=rand, ymin = low.97.5, ymax = high.97.5))+
    geom_errorbarh(aes(y=rand.97.5, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i97.5])), nrow=2, ncol=2)
dev.off()

## Figure A16

data_mun_results <- data_results_canton_mun$data_mun_results

png("Plots/Figure A16.png", width=2500, height=2500, res = 300)
ggarrange(
  ggplot(data_mun_results)+
    geom_point(aes(x=rand, y = rand.2.5))+
    geom_errorbar(aes(x=rand, ymin = low.2.5, ymax = high.2.5))+
    geom_errorbarh(aes(y=rand.2.5, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i2.5])),
  
  ggplot(data_mun_results)+
    geom_point(aes(x=rand, y = rand.25))+
    geom_errorbar(aes(x=rand, ymin = low.25, ymax = high.25))+
    geom_errorbarh(aes(y=rand.25, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i25])),
  
  ggplot(data_mun_results)+
    geom_point(aes(x=rand, y = rand.75))+
    geom_errorbar(aes(x=rand, ymin = low.75, ymax = high.75))+
    geom_errorbarh(aes(y=rand.75, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i75])),
  
  ggplot(data_mun_results)+
    geom_point(aes(x=rand, y = rand.97.5))+
    geom_errorbar(aes(x=rand, ymin = low.97.5, ymax = high.97.5))+
    geom_errorbarh(aes(y=rand.97.5, xmin = low, xmax = high))+
    theme_minimal()+
    labs(x=expression(b[i50]))+
    labs(y=expression(b[i97.5])), nrow=2, ncol=2)
dev.off()

